/*=================================================================================
 * The functions of Flowfield
 *=================================================================================*/

/*---------------------------------------------------------------------------------
 * Standard C++ library headers
 *---------------------------------------------------------------------------------*/

#include <iomanip>

/*---------------------------------------------------------------------------------
 * Local headers
 *---------------------------------------------------------------------------------*/

#include "Flowfield.h"
#include "FileReader.h"
#include "Differ.h"
#include "Global.h"
#include "mpi.h"
#include "omp.h"

/*---------------------------------------------------------------------------------
	Name:		Flowfield::InputRead

	Input(4):	The input file; The mesh nodes and boundary nodes;

	Function:	Read the input file and construct the flow field

	Return:		None
 *---------------------------------------------------------------------------------*/

void Flowfield::InputRead(char *initialization, Array<double> xnode, Array<double> ynode, int bc)
{

	this->xnode = xnode;
	this->ynode = ynode;
	this->bc = bc;

	ni = xnode.GetSize() - 2 * bc;
	nj = ynode.GetSize() - 2 * bc;

	FileReader Input;
	Input.readFile(initialization);

	NS = Input.getIntParameter("NS");
	NR = Input.getIntParameter("NR");

	/* Construct the flow field */
	Construction();

	/* Assign the input value */
	P.Fill(Input.getdoubleParameter("P"));
	T.Fill(Input.getdoubleParameter("T"));
	// Ma.Fill(Input.getdoubleParameter("Mach"));
	V3 = Input.getdoubleParameter("V3");
	cita = Input.getdoubleParameter("angle");

	Pjet = Input.getdoubleParameter("Pjet");
	Tjet = Input.getdoubleParameter("Tjet");
	Ujet = Input.getdoubleParameter("Ujet");
	Vjet = Input.getdoubleParameter("Vjet");
	Ljet = Input.getdoubleParameter("Ljet");

	return;
}

/*---------------------------------------------------------------------------------
	Name:		Flowfield::Construction

	Input(0): 	None

	Function:	Construct the flow field

	Return:		None
 *---------------------------------------------------------------------------------*/

void Flowfield::Construction()
{

	U.Initial(ni + 2 * bc, nj + 2 * bc);
	V.Initial(ni + 2 * bc, nj + 2 * bc);
	P.Initial(ni + 2 * bc, nj + 2 * bc);
	D.Initial(ni + 2 * bc, nj + 2 * bc);
	T.Initial(ni + 2 * bc, nj + 2 * bc);
	C.Initial(ni + 2 * bc, nj + 2 * bc);
	Ma.Initial(ni + 2 * bc, nj + 2 * bc);
	Wav.Initial(ni + 2 * bc, nj + 2 * bc);
	Rgas.Initial(ni + 2 * bc, nj + 2 * bc);
	Cp.Initial(ni + 2 * bc, nj + 2 * bc);
	H.Initial(ni + 2 * bc, nj + 2 * bc);
	E.Initial(ni + 2 * bc, nj + 2 * bc);
	Gamma.Initial(ni + 2 * bc, nj + 2 * bc);

	Mr.Initial(ni + 2 * bc, nj + 2 * bc, NS);

	Mc.Initial(ni + 2 * bc, nj + 2 * bc, NS);
	Mi.Initial(ni + 2 * bc, nj + 2 * bc, NS);
	Yi.Initial(ni + 2 * bc, nj + 2 * bc, NS);
	Di.Initial(ni + 2 * bc, nj + 2 * bc, NS);

	Cpi.Initial(NS);
	Hi.Initial(NS);
	Ei.Initial(NS);

	Uint.Initial(nj);
	Vint.Initial(nj);
	Pint.Initial(nj);
	Dint.Initial(nj);
	Tint.Initial(nj);
	Hint.Initial(nj);
	Eint.Initial(nj);
	Gint.Initial(nj);
	Yint.Initial(nj, NS);

	CS.Initial(ni + 2 * bc, nj + 2 * bc, NS + 3);
	F.Initial(ni + 2 * bc, nj + 2 * bc, NS + 3);
	G.Initial(ni + 2 * bc, nj + 2 * bc, NS + 3);

	Yi_temp0.Initial(NS);
	PLR.Initial(ni + 2 * bc, nj + 2 * bc, 2);
	DLR.Initial(ni + 2 * bc, nj + 2 * bc, 2);
	ULR.Initial(ni + 2 * bc, nj + 2 * bc, 2);
	VLR.Initial(ni + 2 * bc, nj + 2 * bc, 2);
	HLR.Initial(ni + 2 * bc, nj + 2 * bc, 2);
	GLR.Initial(ni + 2 * bc, nj + 2 * bc, 2);
	YLR.Initial(ni + 2 * bc, nj + 2 * bc, 2);
	Yi_temp.Initial(ni + 2 * bc, nj + 2 * bc);
	YL_temp.Initial(ni + 2 * bc, nj + 2 * bc, NS);
	YR_temp.Initial(ni + 2 * bc, nj + 2 * bc, NS);
	YL.Initial(NS);
	YR.Initial(NS);

	Partial_T.Initial(ni + 2 * bc, nj + 2 * bc, NS + 3);
	RHS.Initial(ni + 2 * bc, nj + 2 * bc, NS + 3);

	Yjet.Initial(NS);

	// AM.Initial(NS + 3);

	// UC.Initial(ni + 2 * bc, NS + 3, 3);
	// UUC.Initial(ni + 2 * bc, NS + 3, 3);
	// UCC.Initial(ni + 2 * bc, NS + 3);
	// FX.Initial(ni + 2 * bc, NS + 3);
	// EVR.Initial(ni + 2 * bc, NS + 3, NS + 3);
	// EVL.Initial(ni + 2 * bc, NS + 3, NS + 3);
	// AA.Initial(NS + 3);
	// GG.Initial(ni + 2 * bc, NS + 3, 2);
	// FF.Initial(ni + 2 * bc, NS + 3, 2);
	// FFF.Initial(ni + 2 * bc, NS + 3, 2);
	// FHR.Initial(ni + 2 * bc, NS + 3);
	// FHL.Initial(ni + 2 * bc, NS + 3);
	// LUS.Initial(ni + 2 * bc, NS + 3);
	// YM.Initial(NS);
	// H_old.Initial(ni + 2 * bc);
	// Gamma_old.Initial(ni + 2 * bc);

	return;
}

/*---------------------------------------------------------------------------------
	Name:		Flowfield::FieldInitial

	Input(6):	The component gas constant; The component molecular weight;
				The component mole fraction; The component mass fraction;
				The NASA thermochemical polynomial coefficient

	Function:	Initialize the flow field

	Return:		None
 *---------------------------------------------------------------------------------*/

void Flowfield::FieldInitial(Array<double> Ri, Array<double> Mw, Array<double> Mi_temp, Array<double> Yi_temp, Array<double> Coeff0, Array<double> Coeff1)
{
	dt = 0.0;
	time = 0.0;

	this->Ri = Ri;
	this->Mw = Mw;
	this->Coeff0 = Coeff0;
	this->Coeff1 = Coeff1;

	for (int i = bc; i < ni + bc; i++)
		for (int j = bc; j < nj + bc; j++)
		{
			Wav(i, j) = 1.0 / Fun.sum(3, Yi_temp, Mw);
			Rgas(i, j) = R * 1000 / Wav(i, j);
			D(i, j) = P(i, j) / (Rgas(i, j) * T(i, j));
			// T(i, j) = P(i, j) / D(i, j) / Rgas(i, j);
			for (int k = 0; k < NS; k++)
			{
				Mi(i, j, k) = Mi_temp(k);
				Yi(i, j, k) = Yi_temp(k);
				Di(i, j, k) = Yi(i, j, k) * D(i, j);
				Mc(i, j, k) = Di(i, j, k) / Mw(k) * 1000;
				Cpi(k) = React.GetCpi(T(i, j), Ri(k), k, Coeff0, Coeff1);
				Hi(k) = React.GetHi(T(i, j), Ri(k), k, Coeff0, Coeff1);
				Ei(k) = Hi(k) - Ri(k) * T(i, j);
			}
			Cp(i, j) = Fun.sum(2, Cpi, Yi_temp) * 1000;
			H(i, j) = Fun.sum(2, Yi_temp, Hi) * 1000;
			E(i, j) = Fun.sum(2, Yi_temp, Ei) * 1000;
			Gamma(i, j) = Cp(i, j) / (Cp(i, j) - R * Fun.sum(3, Yi_temp, Mw) * 1000);
			C(i, j) = sqrt(Gamma(i, j) * P(i, j) / D(i, j));
			// U(i, j) = Ma(i, j) * C(i, j) * cos(-cita * atan(1) * 4 / 180);
			// V(i, j) = Ma(i, j) * C(i, j) * sin(-cita * atan(1) * 4 / 180);
			U(i, j) = V3 * cos(-cita * atan(1) * 4 / 180);
			V(i, j) = V3 * sin(-cita * atan(1) * 4 / 180);
			Ma(i, j) = sqrt(pow(U(i, j), 2) + pow(V(i, j), 2)) / C(i, j);

			/* Get the parameters of freestream boundary  */
			if (i == bc)
			{
				Uint(j - bc) = U(i, j);
				Vint(j - bc) = V(i, j);
				Pint(j - bc) = P(i, j);
				Dint(j - bc) = D(i, j);
				Tint(j - bc) = T(i, j);
				Hint(j - bc) = H(i, j);
				Eint(j - bc) = E(i, j);
				Gint(j - bc) = Gamma(i, j);
				for (int k = 0; k < NS; k++)
					Yint(j - bc, k) = Yi(i, j, k);
			}

			for (int k = 0; k < NS; k++)
				CS(i, j, k) = D(i, j) * Yi_temp(k);
			CS(i, j, NS) = D(i, j) * U(i, j);
			CS(i, j, NS + 1) = D(i, j) * V(i, j);
			CS(i, j, NS + 2) = D(i, j) * E(i, j) + 0.5 * (pow(U(i, j), 2) + pow(V(i, j), 2)) * D(i, j);
		}

	for (int k = 0; k < NS; k++)
	{
		Yjet(k) = Yint(0, k);
	}

	/* Get the partial T for IMEX method */
	GetPartial_T();
}

/*---------------------------------------------------------------------------------
	Name:		Flowfield::CFLcondition

	Input(2):	The cfl number; The final time

	Function:	Get the time step

	Return:		None
 *---------------------------------------------------------------------------------*/

void Flowfield::CFLcondition_MPI(int myid,double cfl, double Final_Time)
{
	dx = xnode(1) - xnode(0);
	dy = ynode(1) - ynode(0);

	double dt_conv = 1e10, temp = 0.0;

	for (int i = bc; i < ni + bc; i++)
		for (int j = bc; j < nj + bc; j++)
		{
			temp = min(dx / (abs(U(i, j)) + C(i, j)), dy / (abs(V(i, j)) + C(i, j)));
			if (temp < dt_conv)
				dt_conv = temp;
		}
	dt = cfl * dt_conv;

	if (time + dt > Final_Time)
		dt = Final_Time - time;

	/*calculate minimum dt among processes and broadcast to all processes*/
	double dtc = 0;
	MPI_Reduce(&dt, &dtc, 1, MPI_DOUBLE, MPI_MIN, 0, MPI_COMM_WORLD);
	if (myid == 0) dt = dtc;
	MPI_Bcast(&dt, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
	MPI_Barrier(MPI_COMM_WORLD);

	time += dt;
}

/*---------------------------------------------------------------------------------
	Name:		Flowfield::FieldBoundary_MPI

	Input(1):	the id of the process

	Function:	Set the boundary condition

	Return:		None
 *---------------------------------------------------------------------------------*/

void Flowfield::FieldBoundary_MPI(int id,int MPI_size)
{
	/* Left Boundary ---- Freestream */
	
	for (int j = bc; j < nj+bc ; j++)
		for (int i = 0; i < bc; i++)
		{
			U(i, j) = Uint(j - bc);
			V(i, j) = Vint(j - bc);
			P(i, j) = Pint(j - bc);
			D(i, j) = Dint(j - bc);
			T(i, j) = Tint(j - bc);
			H(i, j) = Hint(j - bc);
			Gamma(i, j) = Gint(j - bc);
			E(i, j) = Eint(j - bc);
			for (int k = 0; k < NS; k++)
				Yi(i, j, k) = Yint(j - bc, k);
		}

	/* Right Boundary ---- Outflow */
	for (int j = bc; j < nj+bc ; j++)
		for (int i = ni + bc; i < ni + 2 * bc; i++)
		{
			U(i, j) = U(2 * (ni + bc) - 1 - i, j);
			V(i, j) = V(2 * (ni + bc) - 1 - i, j);
			P(i, j) = P(2 * (ni + bc) - 1 - i, j);
			T(i, j) = T(2 * (ni + bc) - 1 - i, j);

			for (int k = 0; k < NS; k++)
			{
				Yi(i, j, k) = Yi(2 * (ni + bc) - 1 - i, j, k);
				Yi_temp0(k) = Yi(i, j, k);
				Mc(i, j, k) = Mc(2 * (ni + bc) - 1 - i, j, k);
				Cpi(k) = React.GetCpi(T(i, j), Ri(k), k, Coeff0, Coeff1);
				Hi(k) = React.GetHi(T(i, j), Ri(k), k, Coeff0, Coeff1);
				Ei(k) = Hi(k) - Ri(k) * T(i, j);
			}

			Rgas(i, j) = Fun.sum(2, Yi_temp0, Ri) * 1000;
			D(i, j) = P(i, j) / (T(i, j) * Rgas(i, j));
			Cp(i, j) = Fun.sum(2, Cpi, Yi_temp0) * 1000;
			H(i, j) = Fun.sum(2, Yi_temp0, Hi) * 1000;
			E(i, j) = Fun.sum(2, Yi_temp0, Ei) * 1000;
			Gamma(i, j) = Cp(i, j) / (Cp(i, j) - Rgas(i, j));
		}

	/* Down Boundary ---- Slip boundary */
	if (id ==0){
		for (int j = 0; j < bc; j++)
		for (int i = bc; i < ni+bc; i++)
		{
			if ((xnode(i) >= Ljet) && (xnode(i) <= Ljet+0.002))
			{
				U(i, j) = Ujet;
				V(i, j) = Vjet;
				P(i, j) = Pjet;
				T(i, j) = Tjet;
				for (int k = 0; k < NS; k++)
					Yi(i, j, k) = Yjet(k);
			}
			else
			{
				U(i, j) = U(i, 2 * bc - 1 - j);
				V(i, j) = -V(i, 2 * bc - 1 - j);
				P(i, j) = P(i, 2 * bc - 1 - j);
				T(i, j) = T(i, 2 * bc - 1 - j);
				for (int k = 0; k < NS; k++)
					Yi(i, j, k) = Yi(i, 2 * bc - 1 - j, k);
			}

			for (int k = 0; k < NS; k++)
			{
				Yi_temp0(k) = Yi(i, j, k);
				Mc(i, j, k) = Mc(i, 2 * bc - 1 - j, k);
				Cpi(k) = React.GetCpi(T(i, j), Ri(k), k, Coeff0, Coeff1);
				Hi(k) = React.GetHi(T(i, j), Ri(k), k, Coeff0, Coeff1);
				Ei(k) = Hi(k) - Ri(k) * T(i, j);
			}

			Rgas(i, j) = Fun.sum(2, Yi_temp0, Ri) * 1000;
			D(i, j) = P(i, j) / (T(i, j) * Rgas(i, j));
			Cp(i, j) = Fun.sum(2, Cpi, Yi_temp0) * 1000;
			H(i, j) = Fun.sum(2, Yi_temp0, Hi) * 1000;
			E(i, j) = Fun.sum(2, Yi_temp0, Ei) * 1000;
			Gamma(i, j) = Cp(i, j) / (Cp(i, j) - Rgas(i, j));
		}
	}
	

	/* Up Boundary ---- Outflow */
	if (id == MPI_size-1){
		for (int j = nj + bc; j < nj + 2 * bc; j++)
		for (int i = bc; i < ni + bc; i++)
		{
			U(i, j) = U(i, 2 * (nj + bc) - 1 - j);
			V(i, j) = V(i, 2 * (nj + bc) - 1 - j);
			P(i, j) = P(i, 2 * (nj + bc) - 1 - j);
			D(i, j) = D(i, 2 * (nj + bc) - 1 - j);
			T(i, j) = T(i, 2 * (nj + bc) - 1 - j);
			H(i, j) = H(i, 2 * (nj + bc) - 1 - j);
			Gamma(i, j) = Gamma(i, 2 * (nj + bc) - 1 - j);
			E(i, j) = E(i, 2 * (nj + bc) - 1 - j);
			for (int k = 0; k < NS; k++)
				Yi(i, j, k) = Yi(i, 2 * (nj + bc) - 1 - j, k);
		}
	}
	
}

/*Only send data in the y-direction due to the special domain decomposition*/

void Flowfield::MPIdataInit(int MPI_size)
{
	int grid_num = bc*ni ;
	int types = 8+NS ;
	int volume = grid_num*types ;
	send_data.Initial(volume);
	recv_data.Initial(volume);
}

void Flowfield::MPIBoundary(int id,int MPI_size)
{
	/* Set process tag */
	int tag1 = 0, tag2 = 100, tag3 = 200, tag4 = 300, tag5 = 400, tag6 = 500;
	int idx = 0,dest = 0, src = 0;
	// Send from down to up //
	if (id != MPI_size -1){
		// Package //
		for (int j = -1; j <= 0; j++)
			for (int i = bc; i < bc+ni; i++)
				{
					send_data(idx) = D(i, nj - j);
					send_data(idx + ni * bc) = U(i, nj - j);
					send_data(idx + 2 * ni * bc) = V(i, nj - j);
					send_data(idx + 3 * ni * bc) = P(i, nj - j);
					send_data(idx + 4 * ni * bc) = T(i, nj - j);
					send_data(idx + 5 * ni * bc) = H(i, nj - j);
					send_data(idx + 6 * ni * bc) = Gamma(i, nj - j);
					send_data(idx + 7 * ni * bc) = C(i, nj - j);
					for (int s = 0; s < NS; s++)
						send_data(8 * ni * bc + idx * NS + s) = Yi(i, nj - j,s);
					idx++;
				}
	}
	// Non-blocking communication //
	MPI_Request request3[2];
	MPI_Status status3[2];
	if (id <MPI_size-1) dest = id+1;
	else if(id == MPI_size-1) dest =0 ;
	MPI_Isend(&send_data(0), send_data.GetSize(), MPI_DOUBLE, dest, tag3, MPI_COMM_WORLD, &request3[0]);
	if (id>0) src = id-1;
	else if (id == 0) src = MPI_size-1 ;
	MPI_Irecv(&recv_data(0), recv_data.GetSize(), MPI_DOUBLE, src, tag3, MPI_COMM_WORLD, &request3[1]);
	MPI_Waitall(2, request3, status3);
	// Unpack data //
	idx = 0;
	if (id != 0){
		for (int j = 0; j < bc; j++)
			for (int i = bc; i < ni + bc; i++)
				{
					D(i, bc - j - 1) = recv_data(idx);
					U(i, bc - j - 1) = recv_data(idx +  ni * bc);
					V(i, bc - j - 1) = recv_data(idx + 2 *  ni * bc);
					P(i, bc - j - 1) = recv_data(idx + 3 *  ni * bc);
					T(i, bc - j - 1) = recv_data(idx + 4 *  ni * bc);
					H(i, bc - j - 1) = recv_data(idx + 5 *  ni * bc);
					Gamma(i, bc - j - 1) = recv_data(idx + 6 *  ni * bc);
					C(i, bc - j - 1) = recv_data(idx + 7 *  ni * bc);
					for (int s = 0; s < NS; s++)
						Yi(i, bc - j - 1, s) = recv_data(8 *  ni * bc + idx * NS + s);
					idx++;
				}
	}
	// Send from up to down //
	// Package //
	idx = 0;
	if (id>0){
		for (int j = -1; j < bc-1; j++)
			for (int i = bc; i < bc+ni; i++)
				{
					send_data(idx) = D(i, bc + j +1);
					send_data(idx +  ni * bc) = U(i, bc + j + 1);
					send_data(idx + 2 *  ni * bc) = V(i, bc + j + 1);
					send_data(idx + 3 *  ni * bc) = T(i, bc + j + 1);
					send_data(idx + 4 *  ni * bc) = P(i, bc + j + 1);
					send_data(idx + 5 *  ni * bc) = H(i, bc + j + 1);
					send_data(idx + 6 *  ni * bc) = Gamma(i, bc + j + 1);
					send_data(idx + 7 *  ni * bc) = C(i, bc + j + 1);
					for (int s = 0; s < NS; s++)
						send_data(8 *  ni * bc + idx * NS + s) = Yi(i, bc + j + 1, s);
					idx++;
				}
	}
	// Non-blocking communication //
	MPI_Request request4[2];
	MPI_Status status4[2];
	if (id>0) dest = id-1;
	else dest = MPI_size-1;
	MPI_Isend(&send_data(0), send_data.GetSize(), MPI_DOUBLE, dest, tag4, MPI_COMM_WORLD, &request4[0]);
	if (id<MPI_size-1) src = id+1;
	else src = 0;
	MPI_Irecv(&recv_data(0), recv_data.GetSize(), MPI_DOUBLE, src, tag4, MPI_COMM_WORLD, &request4[1]);
	MPI_Waitall(2, request4, status4);

	// Unpackage //
	idx = 0;
	if (id != MPI_size-1){
		for (int j = 0; j < bc; j++)
			for (int i = bc; i < ni+bc; i++)
				{
					D(i, nj + bc + j) = recv_data(idx);
					U(i, nj + bc + j) = recv_data(idx +  ni * bc);
					V(i, nj + bc + j) = recv_data(idx + 2 *  ni * bc);
					T(i, nj + bc + j) = recv_data(idx + 3 *  ni * bc);
					P(i, nj + bc + j) = recv_data(idx + 4 *  ni * bc);
					H(i, nj + bc + j) = recv_data(idx + 5 *  ni * bc);
					Gamma(i, nj + bc + j) = recv_data(idx + 6 *  ni * bc);
					C(i, nj + bc + j) = recv_data(idx + 7 *  ni * bc);
					for (int s = 0; s < NS; s++)
						Yi(i, nj + bc + j, s) = recv_data(8 *  ni * bc + idx * NS + s);
					idx++;
				}
	}

}



/*---------------------------------------------------------------------------------
	Name:		Flowfield::Update

	Input(2):	The type of time advance scheme;
				The type of the difference scheme;

	Function:	Solve the advection term

	Return:		None
 *---------------------------------------------------------------------------------*/

void Flowfield::Advection(int _TimeAdv, int _Diff)
{
	void (*diff)(int, Array<double>, Array<double>, Array<double> &, Array<double> &, int, int) = &Diff_Initial;

	/* Identify the difference scheme */
	switch (_Diff)
	{
	case 0: /* 1st-order */
		diff = MUSCL_1;
		break;
	case 1: /* 2nd-order */
		diff = MUSCL_2;
		break;
	}

	/* Solve the Flux */
	AUSM(1, F, diff);
	AUSM(2, G, diff);

	/* Identify the time advance scheme */
	switch (_TimeAdv)
	{
	case 0: /* Euler Explicit */
		Time.EE(NS, dt, xnode, ynode, bc, F, G, RHS);
		break;
	case 1: /* 3rd-order TVD Runge-Kutta */
		Time.TVD_RK3(NS, dt, xnode, ynode, bc, F, G, RHS);
		break;
	}
}

/*---------------------------------------------------------------------------------
	Name:		Flowfield::Update_after_Adv

	Input(0):	None

	Function:	Update the flowfield after solving the advection term

	Return:		None
 *---------------------------------------------------------------------------------*/

void Flowfield::Update_after_Adv()
{
	for (int j = bc; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
			for (int k = 0; k < NS + 3; k++)
				CS(i, j, k) = CS(i, j, k) + RHS(i, j, k);

	Update_after_CS();
}

/*---------------------------------------------------------------------------------
	Name:		Flowfield::Explicit

	Input(0):	None

	Function:	Update the flowfield after solving the reaction term and obtaining Di

	Return:		None
 *---------------------------------------------------------------------------------*/

void Flowfield::Explicit()
{
	for (int j = bc; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
		{
			D(i, j) = 0.0;
			for (int k = 0; k < NS; k++)
			{
				CS(i, j, k) = Di(i, j, k);
				D(i, j) += CS(i, j, k);
			}
			CS(i, j, NS) = U(i, j) * D(i, j);
			CS(i, j, NS + 1) = V(i, j) * D(i, j);
			CS(i, j, NS + 2) = D(i, j) * E(i, j) + 0.5 * (pow(U(i, j), 2) + pow(V(i, j), 2)) * D(i, j);
		}

	Update_after_CS();
}

/* Adaptive */
void Flowfield::Explicit(Array<double> Nchem, int step)
{
	for (int j = bc; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
		{
			if (step == 0 || step < Nchem(i, j))
			{
				D(i, j) = 0.0;
				for (int k = 0; k < NS; k++)
				{
					CS(i, j, k) = Di(i, j, k);
					D(i, j) += CS(i, j, k);
				}
				CS(i, j, NS) = U(i, j) * D(i, j);
				CS(i, j, NS + 1) = V(i, j) * D(i, j);
				CS(i, j, NS + 2) = D(i, j) * E(i, j) + 0.5 * (pow(U(i, j), 2) + pow(V(i, j), 2)) * D(i, j);
			}
		}

	Update_after_CS(Nchem, step);
}

/*---------------------------------------------------------------------------------
	Name:		Flowfield::Update_IMEX

	Input(2):	The mass production rate; The diagonal of the Jacobi matrix

	Function:	Update the flowfield after solving the reaction term by IMEX method

	Return:		None
 *---------------------------------------------------------------------------------*/

void Flowfield::Update_IMEX(Array<double> Wi, Array<double> MD)
{
	/* Obtain the conserved variables */
	for (int j = bc; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
		{
			for (int k = 0; k < NS; k++)
				RHS(i, j, k) = RHS(i, j, k) + Wi(i, j, k) * dt;

			for (int k = 0; k < NS + 3; k++)
				CS(i, j, k) = CS(i, j, k) + RHS(i, j, k) / (1.0 - dt * MD(i, j, k, k));
		}

	Update_after_CS();

	/* Get the partial T for IMEX method */
	GetPartial_T();
}

/*---------------------------------------------------------------------------------
	Name:		Flowfield::Update_after_CS

	Input(0):	None

	Function:	Update the flowfield after obtaining the conservation variables

	Return:		None
 *---------------------------------------------------------------------------------*/

void Flowfield::Update_after_CS()
{
	/* Update the parameters after obtaining the conserved variables */
	for (int j = bc; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
		{
			Array<double> Yi_temp1;
			Yi_temp1.Initial(NS);
			D(i, j) = 0;
			for (int k = 0; k < NS; k++)
			{
				D(i, j) += CS(i, j, k);
				Di(i, j, k) = CS(i, j, k);
				Mc(i, j, k) = Di(i, j, k) / Mw(k) * 1000;
			}
			for (int k = 0; k < NS; k++)
			{
				Yi(i, j, k) = CS(i, j, k) / D(i, j);
				Yi_temp1(k) = Yi(i, j, k);
			}
			U(i, j) = CS(i, j, NS) / D(i, j);
			V(i, j) = CS(i, j, NS + 1) / D(i, j);
			E(i, j) = (CS(i, j, NS + 2) - 0.5 * D(i, j) * (pow(U(i, j), 2) + pow(V(i, j), 2))) / D(i, j);
			Wav(i, j) = 1.0 / Fun.sum(3, Yi_temp1, Mw);
			Rgas(i, j) = R * 1000 / Wav(i, j);
			T(i, j) = Get_temp(T(i, j), i, j);
			P(i, j) = T(i, j) * D(i, j) * Rgas(i, j);
			H(i, j) = P(i, j) / D(i, j) + E(i, j);
			for (int k = 0; k < NS; k++)
				Cpi(k) = React.GetCpi(T(i, j), Ri(k), k, Coeff0, Coeff1);
			Cp(i, j) = Fun.sum(2, Cpi, Yi_temp1) * 1000;
			Gamma(i, j) = Cp(i, j) / (Cp(i, j) - R * Fun.sum(3, Yi_temp1, Mw) * 1000);
			C(i, j) = sqrt(Gamma(i, j) * P(i, j) / D(i, j));
			Ma(i, j) = sqrt(pow(U(i, j), 2) + pow(V(i, j), 2)) / C(i, j);
		}
}

/* Adaptive */
void Flowfield::Update_after_CS(Array<double> Nchem, int step)
{
	/* Update the parameters after obtaining the conserved variables */
	for (int j = bc; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
		{
			if (step == 0 || step <= Nchem(i, j))
			{
				D(i, j) = 0;
				for (int k = 0; k < NS; k++)
				{
					D(i, j) += CS(i, j, k);
					Di(i, j, k) = CS(i, j, k);
					Mc(i, j, k) = Di(i, j, k) / Mw(k) * 1000;
				}
				for (int k = 0; k < NS; k++)
				{
					Yi(i, j, k) = CS(i, j, k) / D(i, j);
					Yi_temp0(k) = Yi(i, j, k);
				}
				U(i, j) = CS(i, j, NS) / D(i, j);
				V(i, j) = CS(i, j, NS + 1) / D(i, j);
				E(i, j) = (CS(i, j, NS + 2) - 0.5 * D(i, j) * (pow(U(i, j), 2) + pow(V(i, j), 2))) / D(i, j);
				Wav(i, j) = 1.0 / Fun.sum(3, Yi_temp0, Mw);
				Rgas(i, j) = R * 1000 / Wav(i, j);
				T(i, j) = Get_temp(T(i, j), i, j);
				P(i, j) = T(i, j) * D(i, j) * Rgas(i, j);
				H(i, j) = P(i, j) / D(i, j) + E(i, j);
				for (int k = 0; k < NS; k++)
					Cpi(k) = React.GetCpi(T(i, j), Ri(k), k, Coeff0, Coeff1);
				Cp(i, j) = Fun.sum(2, Cpi, Yi_temp0) * 1000;
				Gamma(i, j) = Cp(i, j) / (Cp(i, j) - R * Fun.sum(3, Yi_temp0, Mw) * 1000);
				C(i, j) = sqrt(Gamma(i, j) * P(i, j) / D(i, j));
				Ma(i, j) = sqrt(pow(U(i, j), 2) + pow(V(i, j), 2)) / C(i, j);
			}
		}
}

/*---------------------------------------------------------------------------------
	Name:		Flowfield::AUSM

	Input(3):	The flux direction; The corresponding flux;
				The function of the difference scheme

	Function:	Solve the flux at each direction

	Return:		The flux
 *---------------------------------------------------------------------------------*/

void Flowfield::AUSM(int direction, Array<double> &Fi, void (*Diff)(int, Array<double>, Array<double>, Array<double> &, Array<double> &, int, int))
{
	double PL = 0.0, PR = 0.0, DL = 0.0, DR = 0.0, UL = 0.0, UR = 0.0, VL = 0.0, VR = 0.0, HL = 0.0, HR = 0.0, GL = 0.0, GR = 0.0;
	double CL = 0.0, CR = 0.0, MaL = 0.0, MaR = 0.0, B1 = 0.0, B2 = 0.0, G1 = 0.0, G2 = 0.0;
	double CI = 0.0, MaI = 0.0, VI = 0.0;

	Diff(direction, xnode, ynode, PLR, P, bc, 0);
	Diff(direction, xnode, ynode, DLR, D, bc, 0);
	Diff(direction, xnode, ynode, ULR, U, bc, 0);
	Diff(direction, xnode, ynode, VLR, V, bc, 0);
	Diff(direction, xnode, ynode, HLR, H, bc, 0);
	Diff(direction, xnode, ynode, GLR, Gamma, bc, 0);

	for (int k = 0; k < NS; k++)
	{
		for (int i = 0; i < ni + 2 * bc; i++)
			for (int j = 0; j < nj + 2 * bc; j++)
				Yi_temp(i, j) = Yi(i, j, k);

		Diff(direction, xnode, ynode, YLR, Yi_temp, bc, 0);

		for (int j = 1; j < nj + bc; j++)
			for (int i = 1; i < ni + bc; i++)
			{
				YL_temp(i, j, k) = YLR(i, j, 0);
				YR_temp(i, j, k) = YLR(i, j, 1);
			}
	}

	for (int j = 1; j < nj + bc; j++)
		for (int i = 1; i < ni + bc; i++)
		{

			PL = PLR(i, j, 0);
			PR = PLR(i, j, 1);
			DL = DLR(i, j, 0);
			DR = DLR(i, j, 1);
			UL = ULR(i, j, 0);
			UR = ULR(i, j, 1);
			VL = VLR(i, j, 0);
			VR = VLR(i, j, 1);
			HL = HLR(i, j, 0);
			HR = HLR(i, j, 1);
			GL = GLR(i, j, 0);
			GR = GLR(i, j, 1);

			for (int k = 0; k < NS; k++)
			{
				YL(k) = YL_temp(i, j, k);
				YR(k) = YR_temp(i, j, k);
			}

			CL = sqrt(GL * PL / DL);
			CR = sqrt(GR * PR / DR);
			CI = 0.5 * (CL + CR);

			if (direction == 1)
			{
				MaL = UL / CI;
				MaR = UR / CI;
			}
			else if (direction == 2)
			{
				MaL = VL / CI;
				MaR = VR / CI;
			}

			if (abs(MaL) < 1.0)
			{
				B1 = 0.25 * pow((MaL + 1.0), 2) + 0.125 * pow((pow(MaL, 2) - 1.0), 2);
				G1 = 0.25 * pow((MaL + 1.0), 2) * (2.0 - MaL) + 3.0 / 16.0 * MaL * pow((pow(MaL, 2) - 1.0), 2);
			}
			else
			{
				B1 = 0.5 * (MaL + abs(MaL));
				G1 = B1 / MaL;
			}
			if (abs(MaR) < 1.0)
			{
				B2 = -0.25 * pow((MaR - 1.0), 2) - 0.125 * pow((pow(MaR, 2) - 1.0), 2);
				G2 = 0.25 * pow((MaR - 1.0), 2) * (2.0 + MaR) - 3.0 / 16.0 * MaR * pow((pow(MaR, 2) - 1.0), 2);
			}
			else
			{
				B2 = 0.5 * (MaR - abs(MaR));
				G2 = B2 / MaR;
			}

			MaI = B1 + B2;
			VI = MaI * CI;

			if (direction == 1)
			{
				if (VI >= 0.0)
				{
					for (int k = 0; k < NS; k++)
						Fi(i, j, k) = VI * DL * YL(k);
					Fi(i, j, NS) = VI * DL * UL + G1 * PL + G2 * PR;
					Fi(i, j, NS + 1) = VI * DL * VL;
					Fi(i, j, NS + 2) = VI * (DL * HL + 0.5 * DL * (pow(UL, 2) + pow(VL, 2)));
				}
				else
				{
					for (int k = 0; k < NS; k++)
						Fi(i, j, k) = VI * DR * YR(k);
					Fi(i, j, NS) = VI * DR * UR + G1 * PL + G2 * PR;
					Fi(i, j, NS + 1) = VI * DR * VR;
					Fi(i, j, NS + 2) = VI * (DR * HR + 0.5 * DR * (pow(UR, 2) + pow(VR, 2)));
				}
			}
			else if (direction == 2)
			{
				if (VI >= 0.0)
				{
					for (int k = 0; k < NS; k++)
						Fi(i, j, k) = VI * DL * YL(k);
					Fi(i, j, NS) = VI * DL * UL;
					Fi(i, j, NS + 1) = VI * DL * VL + G1 * PL + G2 * PR;
					Fi(i, j, NS + 2) = VI * (DL * HL + 0.5 * DL * (pow(UL, 2) + pow(VL, 2)));
				}
				else
				{
					for (int k = 0; k < NS; k++)
						Fi(i, j, k) = VI * DR * YR(k);
					Fi(i, j, NS) = VI * DR * UR;
					Fi(i, j, NS + 1) = VI * DR * VR + G1 * PL + G2 * PR;
					Fi(i, j, NS + 2) = VI * (DR * HR + 0.5 * DR * (pow(UR, 2) + pow(VR, 2)));
				}
			}
		}
}

/*---------------------------------------------------------------------------------
	Name:		Flowfield::Get_temp

	Input(3):	The temperature of the point;
				The coordinate of the point

	Function:	Get the temperature by Newton iteration

	Return:		The temperature after iteration
 *---------------------------------------------------------------------------------*/

double Flowfield::Get_temp(double T, int i, int j)
{
	double T0 = T, T_temp = T;
	double temp1 = 0.0, temp2 = 0.0;
	int count = 0;

	while (count < 10)
	{
		T0 = T_temp;
		for (int k = 0; k < NS; k++)
		{
			Cpi(k) = React.GetCpi(T0, Ri(k), k, Coeff0, Coeff1);
			Hi(k) = React.GetHi(T0, Ri(k), k, Coeff0, Coeff1);
		}
		temp1 = (Fun.sum(2, Yi_temp0, Hi) - E(i, j) * 1e-3) - Fun.sum(2, Yi_temp0, Ri) * T0;
		temp2 = Fun.sum(2, Yi_temp0, Cpi) - Fun.sum(2, Yi_temp0, Ri);
		T_temp = T0 - temp1 / temp2;
		if (abs(T_temp - T0) < 1e-6)
			break;
		else
			count = count + 1;
	}

	return T_temp;
}

/*---------------------------------------------------------------------------------
	Name:		Flowfield::GetPartial_T

	Input(0):	The array of Partial_T

	Function:	Get the Partial_T

	Return:		None
 *---------------------------------------------------------------------------------*/

void Flowfield::GetPartial_T()
{
	double CV = 0.0, Hi = 0.0;

	for (int j = bc; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
		{
			CV = Cp(i, j) - Rgas(i, j);
			for (int k = 0; k < NS; k++)
			{
				Hi = React.GetHi(T(i, j), Ri(k), k, Coeff0, Coeff1);
				Partial_T(i, j, k) = (0.5 * (pow(U(i, j), 2) + pow(V(i, j), 2)) - Hi + Ri(k) * 1e3 * T(i, j)) / (CV * D(i, j));
			}
			Partial_T(i, j, NS + 1) = -U(i, j) / (CV * D(i, j));
			Partial_T(i, j, NS + 2) = 1.0 / (CV * D(i, j));
		}
}

/* WENO-RK3 Test */
void Flowfield::SweepX()
{
	for (int j = bc; j < nj + bc; j++)
		for (int step = 0; step <= 2; step++)
		{
			Boundary_X(j);
			DFWENO(j, step);
			cout << "2" << endl;
			RK3(j, step);
		}
}

void Flowfield::Boundary_X(int j)
{
	/* Left Boundary ---- Freestream */
	for (int i = 0; i < bc; i++)
	{
		U(i, j) = Uint(j - bc);
		V(i, j) = Vint(j - bc);
		P(i, j) = Pint(j - bc);
		D(i, j) = Dint(j - bc);
		T(i, j) = Tint(j - bc);
		H(i, j) = Hint(j - bc);
		Gamma(i, j) = Gint(j - bc);
		E(i, j) = Eint(j - bc);
		for (int k = 0; k < NS; k++)
			Yi(i, j, k) = Yint(j - bc, k);
	}

	// for (int i = 0; i < bc; i++)
	// {
	// 	U(i, j) = -U(2 * bc - 1 - i, j);
	// 	V(i, j) = V(2 * bc - 1 - i, j);
	// 	P(i, j) = P(2 * bc - 1 - i, j);
	// 	D(i, j) = D(2 * bc - 1 - i, j);
	// 	T(i, j) = T(2 * bc - 1 - i, j);
	// 	H(i, j) = H(2 * bc - 1 - i, j);
	// 	Gamma(i, j) = Gamma(2 * bc - 1 - i, j);
	// 	C(i, j) = C(2 * bc - 1 - i, j);
	// 	for (int k = 0; k < NS; k++)
	// 		Yi(i, j, k) = Yi(2 * bc - 1 - i, j, k);
	// }

	/* Right Boundary ---- Outflow */
	for (int i = ni + bc; i < ni + 2 * bc; i++)
	{
		U(i, j) = U(2 * (ni + bc) - 1 - i, j);
		V(i, j) = V(2 * (ni + bc) - 1 - i, j);
		P(i, j) = P(2 * (ni + bc) - 1 - i, j);
		D(i, j) = D(2 * (ni + bc) - 1 - i, j);
		T(i, j) = T(2 * (ni + bc) - 1 - i, j);
		H(i, j) = H(2 * (ni + bc) - 1 - i, j);
		Gamma(i, j) = Gamma(2 * (ni + bc) - 1 - i, j);
		C(i, j) = C(2 * (ni + bc) - 1 - i, j);
		for (int k = 0; k < NS; k++)
			Yi(i, j, k) = Yi(2 * (ni + bc) - 1 - i, j, k);
	}
}

void Flowfield::DFWENO(int j, int step)
{
	double EPWENO = 1.0e-12;
	double DF = 1;
	double q0 = 0.0, q1 = 0.0, q2 = 0.0;
	double D0 = 0.0, D1 = 0.0, D2 = 0.0, D3 = 0.0;
	double IS0 = 0.0, IS1 = 0.0, IS2 = 0.0;
	double a0 = 0.0, a1 = 0.0, a2 = 0.0;
	double W0 = 0.0, W1 = 0.0, W2 = 0.0;
	EigenvalueX(j, step, DF);
	if (P.IsNan())
	{
		cout << "P is nan\n";
		abort();
	}
	cout << "start" << endl;
	for (int i = bc - 1; i < ni + bc - 1; i++)
	{
		if (step == 0)
		{
		}
		for (int n = i - 2; n <= i + 3; n++)
		{
			UC(n, NS + 2, step) = P(n, j) / (Gamma(i + 1, j) - 1.0) + 0.5 * D(n, j) * (pow(U(n, j), 2) + pow(V(n, j), 2));
			// cout << P(n, j) << '\t' << Gamma(i + 1, j) << endl;
			// cout << UC(n, NS + 2, step) << endl;
			FX(n, NS + 2) = (UC(n, NS + 2, step) + P(n, j)) * U(n, j);
		}

		for (int k = 0; k < NS + 3; k++)
			AM(k) = 1.0e-15;
		for (int n = i - 2; n <= i + 3; n++)
		{
			C(n, j) = sqrt(Gamma(i + 1, j) * P(n, j) / D(n, j));
			AM(1) = max(1.0e-15, abs(U(n, j)));
			AM(NS + 1) = max(1.0e-15, abs(U(n, j) + C(n, j)));
			AM(NS + 2) = max(1.0e-15, abs(U(n, j) - C(n, j)));
		}
		for (int k = 0; k < NS + 1; k++)
			AM(k) = AM(1) * 1.3;
		AM(NS + 1) = AM(NS + 1) * 1.3;
		AM(NS + 2) = AM(NS + 2) * 1.3;
		/* Calculate F+ */
		for (int k = 0; k < NS + 3; k++)
			for (int n = i - 2; n <= i + 2; n++)
			{
				for (int m = 0; m < NS + 3; m++)
				{
					AA(m) = (FX(n, m) + AM.MaxValue() * UC(n, m, step)) / 2.0;
					// cout<<AA(m)<<endl;
					GG(n, k, 0) += EVL(i, k, m) * AA(m);
				}
				// cout<<GG(n, k, 0)<<endl;
			}

		for (int k = 0; k < NS + 3; k++)
		{
			q0 = 1.0 / 3.0 * GG(i - 2, k, 0) - 7.0 / 6.0 * GG(i - 1, k, 0) + 11.0 / 6.0 * GG(i, k, 0);
			q1 = -1.0 / 6.0 * GG(i - 1, k, 0) + 5.0 / 6.0 * GG(i, k, 0) + 1.0 / 3.0 * GG(i + 1, k, 0);
			q2 = 1.0 / 3.0 * GG(i, k, 0) + 5.0 / 6.0 * GG(i + 1, k, 0) - 1.0 / 6.0 * GG(i + 2, k, 0);
			D0 = GG(i - 2, k, 0) - GG(i - 1, k, 0);
			D1 = GG(i - 1, k, 0) - GG(i, k, 0);
			D2 = GG(i, k, 0) - GG(i + 1, k, 0);
			D3 = GG(i + 1, k, 0) - GG(i + 2, k, 0);
			IS0 = 13.0 / 12.0 * pow(D0 - D1, 2) + 1.0 / 4.0 * pow(D0 - 3.0 * D1, 2);
			IS1 = 13.0 / 12.0 * pow(D1 - D2, 2) + 1.0 / 4.0 * pow(D1 + D2, 2);
			IS2 = 13.0 / 12.0 * pow(D2 - D3, 2) + 1.0 / 4.0 * pow(3.0 * D2 - D3, 2);
			a0 = 1.0 / 10.0 / pow(EPWENO + IS0, 2);
			a1 = 6.0 / 10.0 / pow(EPWENO + IS1, 2);
			a2 = 3.0 / 10.0 / pow(EPWENO + IS2, 2);
			W0 = a0 / (a0 + a1 + a2);
			W1 = a1 / (a0 + a1 + a2);
			W2 = a2 / (a0 + a1 + a2);
			FF(i, k, 0) = W0 * q0 + W1 * q1 + W2 * q2;
		}

		for (int k = 0; k < NS + 3; k++)
			for (int m = 0; m < NS + 3; m++)
				FFF(i, k, 0) += EVR(i, k, m) * FF(i, m, 0);

		/* Calculate F- */
		for (int k = 0; k < NS + 3; k++)
			for (int n = i - 1; n <= i + 3; n++)
			{
				for (int m = 0; m < NS + 3; m++)
				{
					AA(m) = (FX(n, m) - AM.MaxValue() * UC(n, m, step)) / 2.0;
					GG(n, k, 1) += EVL(i, k, m) * AA(m);
				}
			}

		for (int k = 0; k < NS + 3; k++)
		{
			q0 = -1.0 / 6.0 * GG(i - 1, k, 1) + 5.0 / 6.0 * GG(i, k, 1) + 1.0 / 3.0 * GG(i + 1, k, 1);
			q1 = 1.0 / 3.0 * GG(i, k, 1) + 5.0 / 6.0 * GG(i + 1, k, 1) - 1.0 / 6.0 * GG(i + 2, k, 1);
			q2 = 11.0 / 6.0 * GG(i + 1, k, 1) - 7.0 / 6.0 * GG(i + 2, k, 1) + 1.0 / 3.0 * GG(i + 3, k, 1);
			D0 = GG(i - 1, k, 1) - GG(i, k, 1);
			D1 = GG(i, k, 1) - GG(i + 1, k, 1);
			D2 = GG(i + 1, k, 1) - GG(i + 2, k, 1);
			D3 = GG(i + 2, k, 1) - GG(i + 3, k, 1);
			IS0 = 13.0 / 12.0 * pow(D0 - D1, 2) + 1.0 / 4.0 * pow(D0 - 3.0 * D1, 2);
			IS1 = 13.0 / 12.0 * pow(D1 - D2, 2) + 1.0 / 4.0 * pow(D1 + D2, 2);
			IS2 = 13.0 / 12.0 * pow(D2 - D3, 2) + 1.0 / 4.0 * pow(3.0 * D2 - D3, 2);
			a0 = 3.0 / 10.0 / pow(EPWENO + IS0, 2);
			a1 = 6.0 / 10.0 / pow(EPWENO + IS1, 2);
			a2 = 1.0 / 10.0 / pow(EPWENO + IS2, 2);
			W0 = a0 / (a0 + a1 + a2);
			W1 = a1 / (a0 + a1 + a2);
			W2 = a2 / (a0 + a1 + a2);
			FF(i, k, 1) = W0 * q0 + W1 * q1 + W2 * q2;
		}

		for (int k = 0; k < NS + 3; k++)
			for (int m = 0; m < NS + 3; m++)
				FFF(i, k, 1) += EVR(i, k, m) * FF(i, m, 1);
	}

	for (int i = 2; i < ni + bc - 1; i++)
		for (int k = 0; k < NS + 3; k++)
			FHR(i, k) = FFF(i, k, 0) + FFF(i, k, 1);

	DF = 0;
	EigenvalueX(j, step, DF);
	for (int i = bc; i < ni + bc; i++)
	{
		if (step == 0)
		{
		}
		for (int n = i - 2; n <= i + 3; n++)
		{
			UC(n, NS + 2, step) = P(n, j) / (Gamma(i, j) - 1.0) + 0.5 * D(n, j) * (pow(U(n, j), 2) + pow(V(n, j), 2));
			FX(n, NS + 2) = (UC(n, NS + 2, step) + P(n, j)) * U(n, j);
		}
		for (int k = 0; k < NS + 3; k++)
			AM(k) = 1.0e-15;
		for (int n = i - 2; n <= i + 3; n++)
		{
			C(n, j) = sqrt(Gamma(i, j) * P(n, j) / D(n, j));
			AM(1) = max(1.0e-15, abs(U(n, j)));
			AM(NS + 1) = max(1.0e-15, abs(U(n, j) + C(n, j)));
			AM(NS + 2) = max(1.0e-15, abs(U(n, j) - C(n, j)));
		}
		for (int k = 0; k < NS + 1; k++)
			AM(k) = AM(1) * 1.3;
		AM(NS + 1) = AM(NS + 1) * 1.3;
		AM(NS + 2) = AM(NS + 2) * 1.3;
		/* Calculate F+ */
		for (int k = 0; k < NS + 3; k++)
			for (int n = i - 2; n <= i + 2; n++)
			{
				for (int m = 0; m < NS + 3; m++)
				{
					AA(m) = (FX(n, m) + AM.MaxValue() * UC(n, m, step)) / 2.0;
					GG(n, k, 0) += EVL(i, k, m) * AA(m);
				}
			}

		for (int k = 0; k < NS + 3; k++)
		{
			q0 = 1.0 / 3.0 * GG(i - 2, k, 0) - 7.0 / 6.0 * GG(i - 1, k, 0) + 11.0 / 6.0 * GG(i, k, 0);
			q1 = -1.0 / 6.0 * GG(i - 1, k, 0) + 5.0 / 6.0 * GG(i, k, 0) + 1.0 / 3.0 * GG(i + 1, k, 0);
			q2 = 1.0 / 3.0 * GG(i, k, 0) + 5.0 / 6.0 * GG(i + 1, k, 0) - 1.0 / 6.0 * GG(i + 2, k, 0);
			D0 = GG(i - 2, k, 0) - GG(i - 1, k, 0);
			D1 = GG(i - 1, k, 0) - GG(i, k, 0);
			D2 = GG(i, k, 0) - GG(i + 1, k, 0);
			D3 = GG(i + 1, k, 0) - GG(i + 2, k, 0);
			IS0 = 13.0 / 12.0 * pow(D0 - D1, 2) + 1.0 / 4.0 * pow(D0 - 3.0 * D1, 2);
			IS1 = 13.0 / 12.0 * pow(D1 - D2, 2) + 1.0 / 4.0 * pow(D1 + D2, 2);
			IS2 = 13.0 / 12.0 * pow(D2 - D3, 2) + 1.0 / 4.0 * pow(3.0 * D2 - D3, 2);
			a0 = 1.0 / 10.0 / pow(EPWENO + IS0, 2);
			a1 = 6.0 / 10.0 / pow(EPWENO + IS1, 2);
			a2 = 3.0 / 10.0 / pow(EPWENO + IS2, 2);
			W0 = a0 / (a0 + a1 + a2);
			W1 = a1 / (a0 + a1 + a2);
			W2 = a2 / (a0 + a1 + a2);
			FF(i, k, 0) = W0 * q0 + W1 * q1 + W2 * q2;
		}

		for (int k = 0; k < NS + 3; k++)
			for (int m = 0; m < NS + 3; m++)
				FFF(i, k, 0) += EVR(i, k, m) * FF(i, m, 0);

		/* Calculate F- */
		for (int k = 0; k < NS + 3; k++)
			for (int n = i - 1; n <= i + 3; n++)
			{
				for (int m = 0; m < NS + 3; m++)
				{
					AA(m) = (FX(n, m) - AM.MaxValue() * UC(n, m, step)) / 2.0;
					GG(n, k, 1) += EVL(i, k, m) * AA(m);
				}
			}

		for (int k = 0; k < NS + 3; k++)
		{
			q0 = -1.0 / 6.0 * GG(i - 1, k, 1) + 5.0 / 6.0 * GG(i, k, 1) + 1.0 / 3.0 * GG(i + 1, k, 1);
			q1 = 1.0 / 3.0 * GG(i, k, 1) + 5.0 / 6.0 * GG(i + 1, k, 1) - 1.0 / 6.0 * GG(i + 2, k, 1);
			q2 = 11.0 / 6.0 * GG(i + 1, k, 1) - 7.0 / 6.0 * GG(i + 2, k, 1) + 1.0 / 3.0 * GG(i + 3, k, 1);
			D0 = GG(i - 1, k, 1) - GG(i, k, 1);
			D1 = GG(i, k, 1) - GG(i + 1, k, 1);
			D2 = GG(i + 1, k, 1) - GG(i + 2, k, 1);
			D3 = GG(i + 2, k, 1) - GG(i + 3, k, 1);
			IS0 = 13.0 / 12.0 * pow(D0 - D1, 2) + 1.0 / 4.0 * pow(D0 - 3.0 * D1, 2);
			IS1 = 13.0 / 12.0 * pow(D1 - D2, 2) + 1.0 / 4.0 * pow(D1 + D2, 2);
			IS2 = 13.0 / 12.0 * pow(D2 - D3, 2) + 1.0 / 4.0 * pow(3.0 * D2 - D3, 2);
			a0 = 3.0 / 10.0 / pow(EPWENO + IS0, 2);
			a1 = 6.0 / 10.0 / pow(EPWENO + IS1, 2);
			a2 = 1.0 / 10.0 / pow(EPWENO + IS2, 2);
			W0 = a0 / (a0 + a1 + a2);
			W1 = a1 / (a0 + a1 + a2);
			W2 = a2 / (a0 + a1 + a2);
			FF(i, k, 1) = W0 * q0 + W1 * q1 + W2 * q2;
		}

		for (int k = 0; k < NS + 3; k++)
			for (int m = 0; m < NS + 3; m++)
				FFF(i, k, 1) += EVR(i, k, m) * FF(i, m, 1);
	}

	for (int i = bc; i < ni + bc; i++)
		for (int k = 0; k < NS + 3; k++)
			FHL(i, k) = FFF(i, k, 0) + FFF(i, k, 1);

	for (int i = bc; i < ni + bc; i++)
		for (int k = 0; k < NS + 3; k++)
		{
			LUS(i, k) = FHL(i, k) - FHR(i - 1, k);
			// cout << LUS(i, k);
		}
}

void Flowfield::RK3(int j, int step)
{
	for (int k = 0; k < NS + 3; k++)
		for (int i = bc; i < ni + bc; i++)
		{
			if (step == 0)
			{
				UUC(i, k, 1) = UUC(i, k, 0) - dt * LUS(i, k) / dx;
				UCC(i, k) = UUC(i, k, 1);
			}
			else if (step == 1)
			{
				UUC(i, k, 2) = 0.75 * UUC(i, k, 0) + 0.25 * UUC(i, k, 1) - 0.25 * dt * LUS(i, k) / dx;
				UCC(i, k) = UUC(i, k, 2);
			}
			else
			{
				UUC(i, k, 0) = 1.0 / 3.0 * UUC(i, k, 0) + 2.0 / 3.0 * UUC(i, k, 2) - 2.0 / 3.0 * dt * LUS(i, k) / dx;
				UCC(i, k) = UUC(i, k, 0);
			}
		}

	if (step == 0)
	{
		for (int i = bc; i < ni + bc; i++)
		{
			H_old(i) = D(i, j) * H(i, j) - P(i, j) - 0.5 * D(i, j) * (pow(U(i, j), 2) + pow(V(i, j), 2)) - P(i, j) / (Gamma(i, j) - 1.0);
			Gamma_old(i) = Gamma(i, j);
		}
	}
	for (int i = bc; i < ni + bc; i++)
	{
		D(i, j) = 0.0;
		for (int k = 0; k < NS; k++)
		{
			Di(i, j, k) = UCC(i, k);
			D(i, j) += Di(i, j, k);
			Mc(i, j, k) = Di(i, j, k) / Mw(k) * 1000;
		}
		for (int k = 0; k < NS; k++)
		{
			Yi(i, j, k) = UCC(i, k) / D(i, j);
			Yi_temp0(k) = Yi(i, j, k);
		}
		U(i, j) = UCC(i, NS) / D(i, j);
		V(i, j) = UCC(i, NS + 1) / D(i, j);
		P(i, j) = (Gamma_old(i) - 1.0) * (UCC(i, NS + 2) - H_old(i) - D(i, j) * (pow(U(i, j), 2) + pow(V(i, j), 2)) / 2.0);
		Wav(i, j) = 1.0 / Fun.sum(3, Yi_temp0, Mw);
		Rgas(i, j) = R * 1000 / Wav(i, j);
		T(i, j) = P(i, j) / D(i, j) / Rgas(i, j);
		E(i, j) = UCC(i, NS + 2) / D(i, j) - 0.5 * (pow(U(i, j), 2) + pow(V(i, j), 2));
		// H(i, j) = P(i, j) / D(i, j) + E(i, j);
		for (int k = 0; k < NS; k++)
		{
			Cpi(k) = React.GetCpi(T(i, j), Ri(k), k, Coeff0, Coeff1);
			Hi(k) = React.GetHi(T(i, j), Ri(k), k, Coeff0, Coeff1);
		}
		Cp(i, j) = Fun.sum(2, Cpi, Yi_temp0) * 1000;
		if (step == 2)
		{
			Gamma(i, j) = Cp(i, j) / (Cp(i, j) - R * Fun.sum(3, Yi_temp0, Mw) * 1000);
			H(i, j) = Fun.sum(2, Hi, Yi_temp0) * 1000 + 0.5 * (pow(U(i, j), 2) + pow(V(i, j), 2));
		}
		else
		{
			Gamma(i, j) = Gamma_old(i);
			H(i, j) = H_old(i) / D(i, j) + 0.5 * (pow(U(i, j), 2) + pow(V(i, j), 2)) + P(i, j) / D(i, j) * Gamma(i, j) / (Gamma(i, j) - 1.0);
		}
		C(i, j) = sqrt(Gamma(i, j) * P(i, j) / D(i, j));
		Ma(i, j) = sqrt(pow(U(i, j), 2) + pow(V(i, j), 2)) / C(i, j);
	}
}

void Flowfield::EigenvalueX(int j, int step, int DF)
{
	for (int i = 0; i < ni + 2 * bc; i++)
	{
		for (int k = 0; k < NS; k++)
			UC(i, k, step) = D(i, j) * Yi(i, j, k);
		UC(i, NS, step) = D(i, j) * U(i, j);
		UC(i, NS + 1, step) = D(i, j) * V(i, j);
		UC(i, NS + 2, step) = D(i, j) * H(i, j) - P(i, j);
	}

	if (step == 0)
	{
		for (int i = 0; i < ni + 2 * bc; i++)
		{
			for (int k = 0; k < NS; k++)
				UUC(i, k, step) = D(i, j) * Yi(i, j, k);
			UUC(i, NS, step) = D(i, j) * U(i, j);
			UUC(i, NS + 1, step) = D(i, j) * V(i, j);
			UUC(i, NS + 2, step) = D(i, j) * H(i, j) - P(i, j);
		}
	}

	/*AM(1) = 1.0e-15;
	AM(NS + 2) = 1.0e-15;
	AM(NS + 3) = 1.0e-15;*/

	for (int i = 0; i < ni + 2 * bc; i++)
	{
		for (int k = 0; k < NS; k++)
			FX(i, k) = UC(i, k, step) * U(i, j);
		FX(i, NS) = UC(i, NS, step) * U(i, j) + P(i, j);
		FX(i, NS + 1) = UC(i, NS, step) * V(i, j);
		FX(i, NS + 2) = (UC(i, NS + 2, step) + P(i, j)) * U(i, j);
		// W(I) = SQRT(D(I))
		AM(1) = max(1.0e-15, abs(U(i, j)));
		AM(NS + 1) = max(1.0e-15, abs(U(i, j) + C(i, j)));
		AM(NS + 2) = max(1.0e-15, abs(U(i, j) - C(i, j)));
	}
	for (int k = 0; k < NS + 1; k++)
		AM(k) = AM(1) * 1.3;
	AM(NS + 1) = AM(NS + 1) * 1.3;
	AM(NS + 2) = AM(NS + 2) * 1.3;
	double T0 = 0.0, T1 = 0.0, UXM = 0.0, VYM = 0.0, HM = 0.0, GA = 0.0, QM = 0.0, CM = 0.0;

	for (int i = 0; i < ni + bc; i++)
	{
		T0 = sqrt(D(i, j)) / (sqrt(D(i, j)) + sqrt(D(i + 1, j))); // ROE average
		T1 = 1.0 - T0;
		UXM = T0 * U(i, j) + T1 * U(i + 1, j);
		VYM = T0 * V(i, j) + T1 * V(i + 1, j);
		HM = T0 * H(i, j) + T1 * H(i + 1, j);
		for (int k = 0; k < NS; k++)
			YM(k) = T0 * Yi(i, j, k) + T1 * Yi(i + 1, j, k);
		if (DF == 1)
			GA = Gamma(i + 1, j);
		else
			GA = Gamma(i, j);
		QM = 0.5 * (pow(UXM, 2) + pow(VYM, 2));
		cout << (GA - 1.0) * (HM - QM) << endl;
		CM = sqrt((GA - 1.0) * (HM - QM)); // GAM的平均，YI()的平均
		// cout << CM << endl;
		T0 = UXM * CM;
		for (int k = 0; k < NS; k++)
			for (int m = 0; m < NS; m++)
			{
				if (k == m)
					EVR(i, k, m) = 1.0;
				else
					EVR(i, k, m) = 0.0;
			}
		for (int k = 0; k < NS; k++)
		{
			EVR(i, NS, k) = UXM;
			EVR(i, NS + 1, k) = VYM;
			EVR(i, NS + 2, k) = QM;
			EVR(i, k, NS) = 0.0;
			EVR(i, k, NS + 1) = YM(k);
			EVR(i, k, NS + 2) = YM(k);
		}

		EVR(i, NS, NS) = 0.0;
		EVR(i, NS, NS + 1) = UXM + CM;
		EVR(i, NS, NS + 2) = UXM - CM;
		EVR(i, NS + 1, NS) = 1.0;
		EVR(i, NS + 1, NS + 1) = VYM;
		EVR(i, NS + 1, NS + 2) = VYM;
		EVR(i, NS + 2, NS) = VYM;
		EVR(i, NS + 2, NS + 1) = HM + T0;
		EVR(i, NS + 2, NS + 2) = HM - T0;

		double RCM = 0.0, GAM1 = 0.0, B = 0.0;
		RCM = 1.0 / pow(CM, 2);
		GAM1 = GA - 1.0;
		B = GAM1 * QM;
		for (int k = 0; k < NS; k++)
			for (int m = 0; m < NS; m++)
			{
				EVL(i, k, m) = -RCM * YM(k) * B;
				if (k == m)
					EVL(i, k, m) = RCM * (CM * CM - YM(k) * B);
			}
		for (int k = 0; k < NS; k++)
		{
			EVL(i, NS, k) = -RCM * VYM * pow(CM, 2.0);
			EVL(i, NS + 1, k) = 0.5 * RCM * (B - UXM * CM);
			EVL(i, NS + 2, k) = 0.5 * RCM * (B + UXM * CM);
			EVL(i, k, NS) = RCM * YM(k) * GAM1 * UXM;
			EVL(i, k, NS + 1) = RCM * YM(k) * GAM1 * VYM;
			EVL(i, k, NS + 2) = -RCM * YM(k) * GAM1;
		}
		EVL(i, NS, NS) = 0.0;
		EVL(i, NS, NS + 1) = pow(CM, 2.0) * RCM;
		EVL(i, NS, NS + 2) = 0.0;
		EVL(i, NS + 1, NS) = -0.5 * RCM * (GAM1 * UXM - CM);
		EVL(i, NS + 1, NS + 1) = -0.5 * RCM * GAM1 * VYM;
		EVL(i, NS + 1, NS + 2) = 0.5 * RCM * GAM1;
		EVL(i, NS + 2, NS) = -0.5 * RCM * (GAM1 * UXM + CM);
		EVL(i, NS + 2, NS + 1) = -0.5 * RCM * GAM1 * VYM;
		EVL(i, NS + 2, NS + 2) = 0.5 * RCM * GAM1;
	}
}

void Flowfield::EigenvalueY()
{
	for (int j = 0; j < nj + 2 * bc; j++)
		for (int i = 0; i < ni + 2 * bc; i++)
		{
			for (int k = 0; k < NS; k++)
				F(i, j, k) = D(i, j) * Yi(i, j, k);
			F(i, j, NS) = D(i, j) * U(i, j);
			F(i, j, NS + 1) = D(i, j) * U(i, j);
		}
}

void Flowfield::SweepY()
{
}

void Flowfield::Boundary_Y()
{
}

double Flowfield::MPI_TEST()
{
	double sum = 0.0;
	bc = 2;
	ni = 100 ;
	nj = 100 ;
	for (int j = bc; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
		{
			sum = sum+(i+j)/2.0;
		}
	return sum ;
}